/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

qdd: A Quantum Computer C++ Library constructed using BDD's

Copyright (C) 1999-2007 David Greve <TheBeaNerd@gmail.com>

This library is free software; you can redistribute it and/or modify
it under the terms of the GNU Library General Public License as
published by the Free Software Foundation; either version 2 of the
License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
USA.

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

#include "qdd.h"

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumException

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

QuantumException::QuantumException (const char *string) {

  msg = string;

}

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumCondition

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

QuantumCondition::QuantumCondition () {

  QDD_ECHO ( cout << "QuantumCondition()\n" ) ;

  Value = bddfalse;
  CumulativeIndex = bddfalse;

}

QuantumCondition & QuantumCondition::operator & (const QuantumBit & arg) {

  QDD_ECHO ( cout << "cond = cond & qbit\n" ) ;

  qdd_assert(bdd_simplify(CumulativeIndex,arg.IndexBDD) == bddfalse,
	 "cond & qbit : Overlapping Indicies");

  Value           &= arg.Value();
  CumulativeIndex |= arg.IndexBDD;

  return(* this) ;

}

QuantumCondition & QuantumCondition::operator & (QuantumCondition & arg) {

  QDD_ECHO ( cout << "cond = cond & cond\n" ) ;

  qdd_assert(bdd_simplify(arg.CumulativeIndex,CumulativeIndex) == bddfalse,
	 "cond & cond : Overlapping Indicies");

  Value           &= arg.Value ;
  CumulativeIndex |= arg.CumulativeIndex ;

  arg.Value           = bddfalse;
  arg.CumulativeIndex = bddfalse;

  return(* this);

}


#if 0
QuantumCondition::QuantumCondition (const QuantumBit & arg) {

  QDD_ECHO ( cout << "cond = qbit" ) ;

  Value = arg.Value() ;
  CumulativeIndex = arg.IndexBDD;

}
#endif

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumBit

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

bdd  QuantumBit::State = bddfalse;
int *QuantumBit::IndexHeap;
int  QuantumBit::NextFreeIndex = -1;

int  QuantumBit::INITIALIZED = 0;
int  QuantumBit::MIN_USER_QBIT = -1;
int  QuantumBit::MAX_QBITS = -1;

// Something here ain't right ..
extern Integer Ipow(long, long);

void qdd_init(int nbits, int megs) {

  qdd_assert(! QuantumBit::INITIALIZED, "qdd_init: Redundant QDD initialization");

  int QDD_BDD_NODENUM = (int) (megs * (2000000.0/64.0)) ;

  int QDD_BDD_CACHESIZE = (int) (QDD_BDD_NODENUM / 10.0) ;

  QuantumBit::MAX_QBITS = 2*(5*nbits + 1);

  if (QDD_BDD_NODENUM < (Ipow(2,nbits) * (int) (2000000.0/(2 << 16)))) {

    cerr << "** Warning : RAM to QBIT ratio dangerously small\n" ;

  }

  QuantumBit::MIN_USER_QBIT = ilog2(QuantumBit::MAX_QBITS);

  bdd_init(QDD_BDD_NODENUM,QDD_BDD_CACHESIZE);
  bdd_setvarnum(QuantumBit::MAX_QBITS);
  bdd_autoreorder(BDD_REORDER_WIN2);
  bdd_reorder_verbose(1);

  QuantumBit::IndexHeap = new int[QuantumBit::MAX_QBITS];

  int i;
  for (i=0;i<QuantumBit::MAX_QBITS;i++) {
    QuantumBit::IndexHeap[i] = i;
  }
  QuantumBit::NextFreeIndex = QuantumBit::MIN_USER_QBIT;
  
  for (i=QuantumBit::MIN_USER_QBIT;i<QuantumBit::MAX_QBITS;i++) {
    bdd_addvarblock(bdd_ithvar(i),1);
  }

  srand(7);
  
  QuantumBit::INITIALIZED = 1;

}

void QuantumBit::QuantumBitInit(int value) {

  QDD_ECHO ( cout << "QuantumBitInit(" << value << ")\n" ) ;

  qdd_assert(((value == 0) | (value == 1)), "QuantumBitInit: Initial QBit value must be 0 or 1");

  qdd_assert(INITIALIZED, "QuantumBitInit: Uninitialized Quantum State");

  qdd_assert((NextFreeIndex < MAX_QBITS),"QuantumBitInit: Exhausted Quantum Bits");

  int index = IndexHeap[NextFreeIndex++];

  Index = index;

  IndexBDD = Int2IndexBDD(index,MIN_USER_QBIT);

  Var = bdd_ithvar(index);

  Set(value);

}

QuantumBit::QuantumBit() {

  QDD_ECHO ( cout << "QuantumBit()\n" ) ;

  QuantumBitInit(0);

}

QuantumBit::QuantumBit(int value) {

  QDD_ECHO ( cout << "QuantumBit(" << value << ")\n" ) ;

  QuantumBitInit(value);

}

QuantumBit::~QuantumBit() {

  QDD_ECHO ( cout << "qbit[" << Index << "].~QuantumBit()\n" ) ;

  // Perform Quantum Garbage Collection

  if (! ((Value() == bddtrue) || (Value() == bddfalse))) {
    cerr << "** Warning : ~qbit() : Collecting Mix()'d Bit\n" ;
    Sample();
  }

  putValue(bddfalse);

  qdd_assert(NextFreeIndex > MIN_USER_QBIT,
	 "~qbit(): Excessive DeAllocation (?)");

  qdd_assert(Index >= MIN_USER_QBIT,
	 "~qbit() : Corrupt Quantum Bit");

  IndexHeap[--NextFreeIndex] = Index;

  Index = -1;

  IndexBDD = bddfalse;

  Var = bddfalse;

}

inline void QuantumBit::putValue(bdd arg) {
  State = bdd_ite(IndexBDD,arg,State) ;
}

inline bdd QuantumBit::Value() const {
  return(bdd_relprod(State,IndexBDD,Int2IndexBDD(-1,MIN_USER_QBIT)));
}

void QuantumBit::operator ~ () {

  QDD_ECHO ( cout << "(~ qbit [" << Index << "])\n" ) ;

  putValue(!(Value()));

}

void QuantumBit::operator ^= (QuantumCondition &arg) {

  QDD_ECHO ( cout << "qbit[" << Index << "] ^= (cond)\n" ) ;

  qdd_assert(bdd_simplify(arg.CumulativeIndex,IndexBDD) == bddfalse,
   "qbit ^= cond : Overlapping Indicies (^=)");

  putValue(Value() ^ arg.Value) ;

  arg.CumulativeIndex = bddfalse;
  arg.Value = bddfalse;

}

void QuantumBit::operator ^= (const QuantumBit &arg) {

  QDD_ECHO ( cout << "qbit[" << Index << "] ^= (qbit[" << arg.Index << "])\n" ) ;

  qdd_assert((! (IndexBDD == arg.IndexBDD)), // I guess this could be strength reduced :)
	 "qbit ^= qbit : Overlapping Indicies");

  putValue(Value() ^ arg.Value()) ;

}

QuantumCondition & QuantumBit::operator ! () {

  QDD_ECHO ( cout << "cond = (! qbit[" << Index << "])\n" ) ;

  qdd_assert((Condition.Value == bddfalse) && (Condition.CumulativeIndex == bddfalse),
	 "! qbit : Overused Condition");

  Condition.Value = (! Value()) ;
  Condition.CumulativeIndex = IndexBDD;

  return(Condition) ;

}

QuantumCondition & QuantumBit::operator & (const QuantumBit & arg) {

  QDD_ECHO ( cout << "cond = (qbit[" << Index << "] & qbit[" << arg.Index << "])\n" ) ;

  qdd_assert((! (IndexBDD == arg.IndexBDD)), // I guess this could be strength reduced :)
	 "qbit & qbit : Overlapping Indicies");

  qdd_assert((Condition.Value == bddfalse) && (Condition.CumulativeIndex == bddfalse),
	 "qbit & qbit : Overused Condition");

  Condition.Value = (Value() & arg.Value()) ;
  Condition.CumulativeIndex = (IndexBDD | arg.IndexBDD) ;

  return(Condition) ;

}

QuantumCondition & QuantumBit::operator & (QuantumCondition & arg) {

  QDD_ECHO ( cout << "cond = qbit[" << Index << "] & cond\n" ) ;

  qdd_assert(bdd_simplify(arg.CumulativeIndex,IndexBDD) == bddfalse,
	 "qbit & cond : Overlapping Indicies");

  qdd_assert((Condition.Value == bddfalse) && (Condition.CumulativeIndex == bddfalse), 
	 "qbit & cond : Overused Condition");

  Condition.Value           = (Value() & arg.Value) ;
  Condition.CumulativeIndex = (IndexBDD | arg.CumulativeIndex) ;

  arg.Value = bddfalse;
  arg.CumulativeIndex = bddfalse;

  return(Condition) ;

}

void QuantumBit::Mix() {

  QDD_ECHO ( cout << "qbit[" << Index << "].Mix()\n" ) ;

  bdd subState = bdd_ite(IndexBDD,bddfalse,State);
  
  qdd_assert(bdd_exist(subState,Var) == subState, "qbit.Mix() : State Polluted by MIXd Bit");
  
  qdd_assert(bdd_support(bdd_exist(Value(),Var)) == bddfalse, "qbit.Mix() : Bit Value Polluted by other MIXd Bits");

  putValue(!(Var) ^ Value()) ;

}

void QuantumBit::unMixed() {

  QDD_ECHO ( cout << "qbit[" << Index << "].unMixed()\n" ) ;

  qdd_assert(bdd_support(bdd_exist(Value(),Var)) == bddfalse, "qbit.unMixed() : Bit Value Polluted by other MIXd Bits");

}

int QuantumBit::Sample() {

  QDD_ECHO ( cout << "qbit[" << Index << "].Sample()" ) ;

  bdd curr = Value();
  bdd set  = curr;

  int value;

  while (1) {

    if (set == bddfalse) {
      value = 0;
      State = dag_satisfy((! curr), State);
      break;
    }
    
    if (set == bddtrue) {
      value = 1;
      State = dag_satisfy(curr, State);
      break;
    }

    if (chose()) {
      set = bdd_high(set);
    }
    else {
      set = bdd_low(set);
    }

  }

  QDD_ECHO ( cout << " = " << value << "\n" ) ;

  return(value);

}

void QuantumBit::Set(int value) {

  QDD_ECHO ( cout << "qbit[" << Index << "].Set(" << value << ")\n" ) ;

  qdd_assert((value == 1) || (value == 0), "qbit.Set(int) takes a boolean value");

  bdd curr = Value();

  if ((curr != bddtrue) & (curr != bddfalse)) {
    cerr << "** Warning : qbit.Set() : Setting a Mix()'d bit\n" ;
    Sample();
  }

  if (value) {
    putValue(bddtrue);
  }
  else {
    putValue(bddfalse);
  }

}

/////////////////////////////////////////////////////////////////////////////////////
/************************************************************************************

  class::QuantumRegister

************************************************************************************/
/////////////////////////////////////////////////////////////////////////////////////

QuantumRegister::QuantumRegister(int size) {

  QDD_ECHO ( cout << "QuantumRegister( " << size << " )\n" ) ;

  qdd_assert( size > 0, "qreg(): Register Size must be greater than zero");

  QuBit = new QuantumBit[size];

  Size = size;

}

QuantumRegister::~QuantumRegister() {

  QDD_ECHO ( cout << "~QuantumRegister()\n" ) ;

  delete [] this->QuBit;

}

QuantumBit & QuantumRegister::operator [] (int index) {

  qdd_assert(((index >= 0) && (index < Size)),"qreg[] : Index out of range");

  return(QuBit[index]) ;

}

Integer QuantumRegister::Sample() {

  QDD_ECHO ( cout << "qreg.Sample()\n" ) ;

  Integer value = 0;

  int i;
  for (i=(Size-1);i>=0;i--) {
    value <<= 1;
    value += QuBit[i].Sample();
  }

  return(value);

}

void QuantumRegister::Set(const Integer & value) {

  QDD_ECHO ( cout << "qreg.Set(" << value << ")\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i].Set(BIT(value,i)) ;
  }

}

void QuantumRegister::operator ^= (const QuantumRegister & arg) {

  QDD_ECHO ( cout << "qreg ^= qreg\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i] ^= arg.QuBit[i];
  }

}

void QuantumRegister::Mix() {

  QDD_ECHO ( cout << "qreg.Mix()\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i].Mix();
  }

}

void QuantumRegister::unMixed() {

  QDD_ECHO ( cout << "qreg.unMixed()\n" ) ;

  int i;
  for (i=0;i<Size;i++) {
    QuBit[i].unMixed();
  }

}

Integer QuantumRegister::Period() {

  QDD_ECHO ( cout << "qreg.Period()\n" ) ;

  bdd s0 = bddtrue;
  Integer v0 = 0;
  int i;
  for (i=(Size-1);i>=0;i--) {
    v0 <<= 1;
    if (bdd_simplify(QuBit[i].Value(),s0) == bddtrue) {
      v0 += 1;
      s0 &= QuBit[i].Value(); // I don't think we need this.
    }
    else {
      s0 &= (! QuBit[i].Value());
    }
  }

  bdd s1 = (! s0);
  Integer v1 = 0;
  for (i=(Size-1);i>=0;i--) {
    v1 <<= 1;
    if (bdd_simplify(QuBit[i].Value(),s1) == bddtrue) {
      v1 += 1;
      s1 &= QuBit[i].Value();
    }
    else {
      s1 &= (! QuBit[i].Value());
    }
  }  

  (* this).Sample();

  Integer period = v1 - v0;

  if (s1 == bddfalse) {
    period = Ipow(2,Size);
  }

  return(period);

}
